function collocInfo = collocateInfo( segs, nodes )
%COLLOCATEINFO Summary of this function goes here
%   Detailed explanation goes here
nseg = size(segs) - 1; % time intervals, default 'segs' = [-1; 1]
nnode_all = sum(nodes);
end_idx = 0;
tau_all = zeros(nnode_all, 1); w_all = zeros(nnode_all, 1);
D_all = zeros(nnode_all, nnode_all + 1); A_all = zeros(nnode_all, nnode_all);
for iseg = 1 : nseg
    [tau, w] = getLGR(nodes(iseg));
    coef0 = (segs(iseg) + segs(iseg+1))/2; coef1 = (segs(iseg+1) - segs(iseg))/2;
    tau_seg = coef0 + tau*coef1;
    D_seg = getDiffMat([tau_seg; segs(iseg+1)]);
    A_seg = inv(D_seg(:,2:end));
    start_idx = end_idx + 1; end_idx = end_idx + nodes(iseg);
    tau_all(start_idx : end_idx) = tau_seg;
    w_all(start_idx : end_idx) = w * coef1;
    D_all(start_idx : end_idx, start_idx : end_idx+1) = D_seg;
    A_all(start_idx : end_idx, start_idx : end_idx) = A_seg;
end
% TODO: form a composite diff matrix for all segments in one phase

% set struct
collocInfo.points = tau_all;
collocInfo.weights = w_all;
collocInfo.diffMat = sparse(D_all);
%collocInfo.diffMatDiag = diag(D_all);

% collocInfo.intMat = sparse(A_all);   % integration matrix

    function [tau, w] = getLGR( nnodes )
        N = nnodes - 1; N1 = N + 1;
        tau = -cos(2*pi*(0:N)/(2*N + 1))'; tau_old = 2; % Initial guess for LGR nodes
        P = zeros(N1,N1 + 1);     % The Legendre Vandermonde Matrix      
        free = 2 : N1;            % Free abscissae
        while max(abs(tau-tau_old)) > 2.22e-15 % 10 * eps(_in_matlab)
            tau_old = tau;
            P(1,:) = (-1).^(0 : N1); P(free,1) = 1; P(free,2) = tau(free);
            for k = 2 : N1, P(free, k+1) = ( (2*k - 1)*tau(free).*P(free, k)-(k-1)*P(free, k-1) )/k; end
            tau(free) = tau_old(free)-((1-tau_old(free))/N1).*(P(free,N1)+P(free,N1+1))./(P(free,N1)-P(free,N1+1));
        end
        P = P(1:N1, 1:N1); % The Legendre-Gauss-Radau Vandermonde
        w = zeros(N1, 1); w(1) = 2/N1^2; % Compute the weights
        w(free) = (1-tau(free))./(N1 * P(free, N1)).^2;
    end
    function D = getDiffMat(tau_aug) %---- Compute collocation differentiation matrix
        M = length(tau_aug); M1 = M+1; M2 = M*M;
        Y = repmat(tau_aug,1,M); Ydiff = Y - Y' + eye(M); % Compute the barycentric weights
        WW = repmat(1./prod(Ydiff,2), 1,M); % product each row
        D = WW./(WW'.*Ydiff);
        D(1:M1:M2) = 1-sum(D); D = -D';     % Diff matrix
        D = D(1:end-1,:);                   % Augment for LGR points
    end
end

